Invasion Percolation with Temperature and the Nature of SOC in Real Systems 
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We show that the introduction of thermal noise in Invasion Percolation (IP) brings the system 
outside the critical point. This result suggests a possible definition of SOC systems as ordinary 
critical systems where the critical point correspond to set to one of the parameters. We recover 
both IP and EDEN model, for T ^ 0, and T ^ oo respectively. For small T we find a dynamical 
second order transition with correlation length diverging when T -^ 0. 
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The spontaneous development of complex and fractal 
structures has been studied on the basis of several mod- 
els which manifest the property of Self-Organized Crit- 
icality (SOC) ||l|. This concept is very intriguing and 
its very meaning has been highly debated. The combi- 
nation of different properties as for example stochastic 
and quenched disorder, usually brings out of the criti- 
cality. Henceforth, the distinction with ordinary critical 
phenomena (instead of self-organized) seems to become 
elusive. In order to clarify these basic questions we con- 
sider here one of the classical models of self-organization, 
the Invasion Percolation (IP) model 0. IP describes the 
displacement of a fluid in a disordered net of random 
throats due to another immiscible fluid pushed with a 
vanishing pressure rate. 

In this letter we study this model when a temperature- 
like noise T is present. The reason for this generalization 
is twofold. On one side, one is interested in the robustness 
of SOC with respect to external solicitations pM| . In par- 
ticular we find that for T ^ a finite correlation length 
appears and the process goes out of criticality. This result 
(togheter with the observation that in other SOC model 
the scaling properties are limited by the "slow driving" 
H) suggests a possible definition for SOC phenomena in 
real systems. A system or a dynamical process is SOC 
if the critical value of the driving parameter is 0, instead 
of another real number. The reason why such a value 
makes such a large difference is because the driving pa- 
rameter of these process is always a ratio (grain of sand 
added with respect to the total number of sites for the 
sandpiles, sites whose "value" is changed for IP, DLA ||] 
Bak and Sneppen 0] etc.) and any value smaller than 
a certain threshold can be considered equal to 0. For 
this reason the zero value tends to occupy a much larger 
region of the phase space than the other real numbers. 

On the other side IP model is the most famous and 
simple example of evolution in quenched disorder. IP 
dynamics of evolution is deterministic and extremal in 
the sense that at each time step the fluid invades the 
environment by selecting the minimum throat available. 
This limiting case is particularly instructive since the ex- 
tremal dynamics is suitable to be studied analytically. By 
using the tool of the Run Time Statistics (RTS) ||,|n|, 
we can approach also the more real stochastic case, where 
fluctuations affect the dynamics of invasion. 



It is useful at this point to describe in some detail the 
IP model that we are studying 0: (a) in a lattice of 
size L a random number Xi (extracted from the uniform 
distribution po{x) = 1 for a; G [0, 1]) is assigned to each 
bond i; (b) at time t = the dynamics starts from a finite 
connected set of bonds Co (in the asymptotic regime the 
system does not depend upon Co). Ct is the connected 
set of bonds grown until time t. At each time step the 
interface dCt of the growing cluster Ct is the set of non- 
grown bonds in contact with the cluster. Only bonds 
belonging to dCt grow at the time step t; (c) At time t 
the bond i £ dCt with the lowest Xi grows; then Ct+i = 
CtU{i} and the interface is updated. The dynamics stops 
when the cluster percolates the lattice. 

This simple growth model develops spontaneously ge- 
ometrical and dynamical critical features. In particular: 
(1) the asymptotic cluster is a fractal (i.e. it has an infi- 
nite correlation length) with fractal dimension Df ~ 1.89 
in a 2d lattice, which is the same fractal dimension of the 
infinite cluster of percolation at the critical point; (2) the 
normalized histogram (jitix), of the interface variables, 
has the following asymptotic shape: 
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while the initial shape is obviously (/'o(a;) = 1. Pc coin- 
cides with the critical threshold of the percolation in the 
lattice; (3) the asymptotic dynamics evolves for critical 
avalanches. Any bond i growing at time t is the initiator 
of an own avalanche. An avalanche is a temporal consec- 
utive sequence of causally and geometrically connected 
growth events starting with the growth of the initiator 
(for a detailed definition of avalanche see e.g. O). Note 
that X of the initiator, due to the shape of the asymptotic 
histogram and of the acceptance profile must he x < Pc 
if t is very large. The size distribution D{s;x) (where x 
is the random number of the initiator) of the avalanche 
has the following behavior: 
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where f(x) — constant for a; ^ 1 and decay exponen- 
tially for a; ^ 1 (i.e. for s > sq ^ \x — pd "'^ ), with 
T = 1.57 ± 0.03 and cr = 1 - T -H 2/Df = 0.49 ± 0.03. 
It may be observed that if x = pc the size distribution 



is a power law and the characteristic size diverges. As 
a consequence, this kind of avalanches are called critical 
avalanches. 

We now generalize the IP model by introducing the 
presence of thermal noise. The numerical study of an 
analogous application to the Bak-Sneppen model can 
be found in pl whilst the case of sandpiles has been con- 
sidered in y7 The first effect of a finite temperature T 
is that the deterministic dynamics becomes stochastic, 
in such a way that the larger the temperature the larger 
the stochasticity. The definitions of Ct and dCt in this 
model are the same of IP, but the growth rule is different: 
each bond i e dCt has the following growth probability 
depending on the realization of the quenched disorder: 
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where (3 — \/T and {a;}act the realization of the 
quenched disorder on the interface dCt- The larger is 
T the more rji^t is independent on Xi. Hereafter, we shall 
indicate with ||Ct|| the number of bonds belonging to Cj, 
and with ||9Ct|| the number of bonds belonging to dCt- 
It is important to study the two different limits T ^ oo 
and T — > 0. In the first limit we have: 
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where l|9Ct|| is the total number of bonds belonging to 
the growth interface at time t. Eq.(Q) means that all the 
bonds on the interface have the same probability to grow. 
This model is well known and usually called Eden model 
l^: this dynamical growth generates a compact cluster 
(fractal dimension equal to the space dimension) with a 
rough surface (interface). In the second limit we have: 
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where dCt — {i} means the interface dCt minus the bond 
i. Eq.(0) provides nothing else the deterministic extremal 
growth rule of IP: rji^t — 1 if a^i is the extremal (minimum) 
value and zero otherwise. In this paper we address mainly 
the study of the behavior for small values of T, i.e. the 
transition of the model towards IP. In particular we will 
study the case of a 2d square bond lattice. We started 
by studying some Monte-Carlo simulations of this model. 
The presence of the temperature introduces a character- 
istic length £,{T) the effect of which is quite clear in Fig.|l| 
where percolating clusters for different values of /? = 1/T 
are shown. The differences between the clusters can be 
explained by characterizing qualitatively the dynamical 
evolution of the growth. 

For any value of T, a characteristic time t*{T) exists 
such that, for t < t*{T) the dynamics of the model is the 
IP dynamics, i.e. even if the dynamical rule given by Eq.|3| 
is not deterministic, the effect of stochasticity is still neg- 
ligible and the effective dynamics is almost extremal. On 



the other hand for t > t* (T) the effect of the stochastic 
noise begins to be more and more important and the de- 
viation from IP and then from fractality, becomes larger. 
If we suppose that T <g; 1, and then t*{T) ^ 1, it is clear 
that t* (T) represents the correlation time of the system. 
Since one bond is removed for each time-step, t*{T) rep- 
resents also the number of bonds so(T') in a correlated re- 
gion of the cluster when t ^ t*{T). This is in agreement 
with the idea that at T > IP is the repulsive fixed point 
of the dynamics under a spatio-temporal coarse-graining 
transformation, and the Eden model is the trivial attrac- 
tive fixed point characterized by T -^ oo. These features 
can be checked by looking at the dynamical evolution of 
the histogram 4>t{x). Obviously (j)o{x) = 1; for i < t*{T) 
as previously noted, the evolution is the same of IP, that 
is (ptix) evolves in the step-function given by Eq.(^. At 
t — t*{T), (j)t{x) is a smoothened step function (the size 
of the smoothened interval around Pc increases with T). 
For t > t*{T), because of stochasticity, the growth of 
bonds with x well larger pc are permitted and the his- 
togram (ptix) shifts towards high values of x. We have 
measured through simulations t*(T) by measuring the 
time step when (ptix) start to shifts and we obtain the 
scaling law t*{T) = sq{T) - T^^ with 7 = 1.9 ± 0.2. In 
the following we find the same behavior theoretically and 
we link it to the correlation length of the structure. 

In order to study analytically the model, we formu- 
late the generalization to stochastic growth dynamics of 
the Run Time Statistics (RTS) [|jl^] that we caU Gen- 
eralized Run Time Statistics (CRTS). The usual RTS is 
a probabilistic technics based on the concept of condi- 
tional probability, introduced to study IP-like dynamics, 
i.e. deterministic extremal dynamics with quenched dis- 
order. With the CRTS approach one can solve the fol- 
lowing problem: suppose to fix the time-ordered path Ct 
followed by the dynamics, and to ignore the realization 
of the disorder: then one can compute the joint prob- 
ability density function Pt{{x}dCt) of ^U the variables 
Xi of the bonds i belonging to the interface dCt, con- 
ditioned to the history Ct- Furthermore one can com- 
pute the conditioned probability of any possible next 
growth step. This joint Probability Density Function 
(PDF) Pt{{x}oct) pla-ys a central role, since from it one 
can compute the probability (conditioned to the whole 
past history, i.e. to all the previous steps of the path) 
of any possible next growth step. After that, one up- 
dates consequently the joint probability density obtain- 
ing Pt+i{{x}gct+i)- Here we expose an approximated 
version of CRTS. The approximation consists in assum- 
ing that at any time-step the PDF can be written as the 
product of single bond density functions Pk,t{xk) 
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This means that one is assuming that all the information 
about the history can be contained in the set of effective 
single bond density functions. Usually this is not the 



case, in fact, even if we start the dynamics with indepen- 
dent variables (as in this case), the information about 
the dynamical history generates correlations among the 
interface variables [|2|. However, it can be seen that this 
approximation works very well even for IP where the the 
effect of this correlation is the maximum, due to the ex- 



tremal nature of the dynamics |13|. 

Starting from the PDF's we want to compute the con- 
ditional probability /ii_t that a certain bond i £ dCt 
grows at time t. Let us suppose to have the "effective" 
probabilities density functions Pk,t{xk) for each k € dCt- 
The functions Pk,t{xk) are determined by the whole past 
history up to time t (in particular, for i = each 
Pk,t{x) = Po{x) = 1 because there is no information yet 
on the dynamics). Starting from functions Pk,t{xk) we 
write the conditioned probability ^i^t as: 
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Eq-(p[) provides the Growth Probability Distribution 
(GPD) conditioned to the past growth history up to time 
t. At this point we may iterate the procedure by updating 
the PDF's, that is by updating the "effective" probability 
density functions conditioned to the last growth-event. In 
order to do that, we have to distinguish three cases: (a) 
the last bond grown i, (b) the other bonds j belonging 
to 9Cf , and finally (c) the bonds just entered in the new 
interface dCt+i because of the growth of i: 
(a) in this case, i does not belong to dCt+i] we introduce 
the new symbol rui^t+iix) analogous to pk,t{x)'. 
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(b) In this case we have: 

Pj,t+i{x) ^ / ••■ / n id.XkPk,t{xk)]- 
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(c) Finally, we have pj^t-^-i{x) — po{x) — 1; let us call 
rii^t the number of these bonds. Note that the following 
relations hold: ||Ct|| = t and ||5Ct+i|| = \\dCt\\ +ni^t - 1- 
Here after we shall call fit and nt the average values over 
different histories respectively of ||Ct|| and rii^f 

Using Eqs.(^,(0),(||) and the rule that bonds just en- 
tered the interface have simply po{x) = 1 as "effective" 
density function, one can describe from a conditional 
probability point of view any possible dynamical his- 
tory, knowing only poix) and the dynamical rule given 
by Eq.(|). In ||l^,|l|| these approach in the T = limit 
has been used to study IP to evaluate both Df and r. 



Now we use this generalized approach to study the tran- 
sition towards IP (stochastic-extremal transition). We 
introduce, then, the histogram ht{x). ht{x) is the distri- 
bution of a;'s on the interface at time t 

ht{x) dx = jj= of interface bonds with throat £ [x,x -[- dx\ 
If wc fix a history up to time t, we can write: 
ht{x) = ^ PiA^) 

iedCt 

where the functions pi^x) must be evaluated through 
the "alghorythm" provided by Eqs.gj^H for the given his- 
tory. Note that J^ dxht{x)—\\dCt\\. Since the disorder is 
quenched (i.e. time-independent), the dynamical equa- 
tion for ht{x) is 



ht+i{x) = ht{x) - mi^t+i{x) + ni^tPo{x) 



(9) 



It is convenient to study the normalized histogram (j)t (x) 
defined as (/)t(a;)=/it(a;)/ II 9Ct II (i.e. /p dx(j)t{x)=l). Since 
(as for IP) (ptix) is an almost self-averaging quantity for 
small T, we can take the average of Eq.|| over all the 
possible histories in order to evaluate it. After some al- 
gebra and some approximations, it is possible to derive 
the following equations: 



nt+if^t+iix) ^^tf^tix) -i^tf^tix) 



1 



l-hr2te''(^-i/«* 



■+nt (10) 



where Ot+i = fit + nt — 1. To obtain Eq.O we have 
assumed that fit ^ I and e^ ^ fit- Clearly the dy- 
namical evolution of the histogram is strictly related to 
that of nt] in IP for t ^ 1 we have nt ~ 1/pc H- For 
t > t * (T) the evolution of (t)t{x) will be very slow (i.e. 
\<j)t+i{x) — (j)t{x)\/<j)t{x) ^ 1), because of the increasing 
effect of stochasticity. The main effect is geometrical and 
it is related to the fact that, as the system pass from frac- 
tal to homogeneous , n^ — > 1 (it can be shown that nt — 1 
represents the asymptotic value of the ratio between the 
number of bonds belonging to the interface and the bonds 
belonging to the cluster). Here we have: 
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(f>t{x) is a smoothened step function centered at x = l/nt 
and large Ax ~ T. For t = t*{T) » 1 we have nt ~ 1/pc, 
as the dynamics is IP-like ]l^,|l^. Then 
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This functions differs from eq. (|l|) only in a region of ex- 
tension Ax ~ T just around x = pc- The agreement 
between this function and the numerical data is very 
good for a wide range of T (Fig.||). From Eq.([l^) and 
from the exponent of IP, we obtain the behavior of so(T) 



and ^(T) for small values of T. In IP an avalanche with 
initiator different from pc of a quantity Ax has a typ- 
ical size so{Ax) ^ Ax~^^'^ . Here we have a natural 
value Ax ~ T even for the maximal sequence of corre- 
lated growth events. Hence so(T) - T"!/'^ = T"'^ with 
7 = 2.0±0.1 in agreement with the simulations. Finally 
because of the fractahty of IP, we have so{T) ^ f (T)^, 
hence CiT) - T""^ with ly = -f/D = 1.10 ± 0.05. 

In conclusion, we presented here a generalization of the 
IP model where stochasticity, by means of a temperature- 
like parameter T is introduced. The model produces 
structures that are fractal and self-organized only by tun- 
ing this parameter to 0, otherwise a finite correlation 
length exists. This behaviour (similar to that observed 
for the Bak and Sneppen model 0] by M. Vergeles [^), 
supports the hypothesis that SOC models are closely re- 
lated to ordinary critical systems, where parameters have 
to be tuned to their critical value. The fundamental dif- 
ference, in our opinion, is in the feasibility of this tuning. 
For SOC models, one has typically to consider limits to 
instead to some other real number, the larger probability 
to achieve this value with respect to any other value is 
linked to the nature of the driving parameter that it is 
usually a density for such systems. 
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FIG. 1. Different percolating clusters for different values of 
/3 = tjt; note that the "fractafity" increases with l3. 
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FiG. 2. Different histograms at i = t*{T) for different 
(3 — y\ the larger is (3 the more IP-like the histogram is. 
Dashed line represents Eq.(|l2[); numerical data are shown by 
fiUed dots. 
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